#*************         NHANES       *************
#*************      Sample Code 01  *************
#*************  使用R代码进行读取   *************
#*************                      *************

#### 示例 ####

#### 0.准备好环境 ####
# Paper01-美国成年人膳食中类胡萝素与认知功能的关系, NHANES 2011-2014
# Paper01-Dietary carotenoids and cognitive function among US adults, NHANES 2011–2014

#下载读取xpt文件的R包“haven”并载入
#install.packages("haven") # 下载一次后无需再次下载，可以在代码前加入 # 将代码转换为注释
library(haven)

#下载用于数据处理包 dplyr, plyr
#install.packages('plyr') # 下载一次后无需再次下载，可以在代码前加入 # 将代码转换为注释
library(plyr)
#install.packages('dplyr') # 下载一次后无需再次下载，可以在代码前加入 # 将代码转换为注释
library(dplyr) # 链接：https://dplyr.tidyverse.org/reference/mutate.html

#下载用于数据快速预览的包 arsenal，使用 tableby 函数
# install.packages('arsenal') # 下载一次后无需再次下载，可以在代码前加入 # 将代码转换为注释
library(arsenal) # 链接：https://cran.r-project.org/web/packages/arsenal/vignettes/tableby.html

#用于加权情况下的分析
library(survey)

#setwd("J:/nhanes/数据分析/NHANES_20221011") #需要转换为自己的数据读取路径
#列出下属的 NHANES 数据文件-list.files()

#### 1. 定位数据模块和变量，获取源数据 ####
##### 1.1 DEMO-人口学数据提取 #####
### 1.1 提取 Component文件
# NHANES官网链接：https://wwwn.cdc.gov/Nchs/Nhanes/2011-2012/DEMO_G.htm
demo.g <- read_xpt("2011-2012/Demographics/demo_g.xpt")#参见上述设置默认路径，文件名称后缀不同"g,h"
demo.h <- read_xpt("2013-2014/Demographics/demo_h.xpt")#参见上述设置默认路径，文件名称后缀不同"g,h"

# demo.h <- read_xpt("J:/nhanes/数据分析/NHANES_20221011/2013-2014/Demographics/demo_h.xpt")#和28行做对比

colnames(demo.g)
colnames(demo.h)

# tableby 函数使用的注意事项-1，分类变量可以在前面加上 factor；如果需要分组，则
tab1 <- tableby( ~ RIDAGEYR + RIAGENDR + RIDRETH3 + DMDEDUC2 + INDFMPIR, 
                 data=demo.g)  #注意！这里是英文的波浪号 ~
summary(tab1, text=TRUE)

tab1 <- tableby( ~ RIDAGEYR + factor(RIAGENDR) + RIDRETH3 + DMDEDUC2 + INDFMPIR, 
                 data=demo.g)   #注意！这里是英文的波浪号 ~
summary(tab1, text=TRUE)

### 1.2 提取研究所需要的变量
# 年龄-RIDAGEYR; 性别-RIAGENDR; 种族-RIDRETH3; 教育程度-DMDEDUC2; 贫困程度-INDFMPIR;

demo.data.file <- dplyr::bind_rows(list(demo.g, demo.h))
dim(demo.data.file)
demo.data <- demo.data.file[,c('SEQN', 'RIDAGEYR', 'RIAGENDR', 'RIDRETH3', 'DMDEDUC2', 'INDFMPIR')]
View(demo.data)#查看数据

tab1 <- tableby( ~ RIDAGEYR + factor(RIAGENDR) + factor(RIDRETH3) + factor(DMDEDUC2) + INDFMPIR, 
                 data=demo.g)  #注意！这里是英文的波浪号 ~
summary(tab1, text=TRUE)
# 
##### 1.2 SMQ-吸烟数据提取 #####
### 1.1 提取 Component文件

smq.g <- read_xpt("2011-2012/Questionnaire/smq_g.xpt")
smq.h <- read_xpt("2013-2014/Questionnaire/smq_h.xpt")
# NHANES官网链接：https://wwwn.cdc.gov/Nchs/Nhanes/2013-2014/SMQ_H.htm

### 2. 提取研究所需要的变量
# 是否吸烟至少100支-SMQ020; 现在是否吸烟-SMQ040; 多久前开始戒烟-SMQ050Q;
# 多久前开始戒烟的时间单位（天、周、月、年）-SMQ050U;

smq.data.file <- dplyr::bind_rows(list(smq.g, smq.h))
colnames(smq.data.file)
smq.data <- smq.data.file[,c('SEQN', 'SMQ020', 'SMQ040', 'SMQ050Q', 'SMQ050U')]

View(smq.data)

##### 1.3 ALQ-饮酒数据提取 ######
### 1.提取 Component文件
# NHANES官网链接：https://wwwn.cdc.gov/Nchs/Nhanes/2013-2014/ALQ_H.htm

alq.g <- read_xpt("2011-2012/Questionnaire/alq_g.xpt")
alq.h <- read_xpt("2013-2014/Questionnaire/alq_h.xpt")

### 2.提取研究所需要的变量
# 每年至少喝12杯酒-ALQ101
# 一生中至少喝过12次酒-ALQ110；
# 过去12个月多久喝一次酒-ALQ120Q
# 过去12个月饮酒频率的单位（周，月，年）-ALQ120U;

alq.data.file <- dplyr::bind_rows(list(alq.g, alq.h))
colnames(alq.data.file)
alq.data <- alq.data.file[,c('SEQN', 'ALQ101', 'ALQ110', 'ALQ120Q', 'ALQ120U')]
View(alq.data)


tab1 <- tableby(~factor(ALQ101) + factor(ALQ110) + ALQ120Q + factor(ALQ120U), 
                data = alq.data)
summary(tab1)

tab1 <- tableby(~factor(ALQ101) + factor(ALQ110) + ALQ120Q + factor(ALQ120U), 
                data = alq.h)
summary(tab1)

##### 1.4 BMX-BMI、腰围数据提取 #####
### 1.提取 Component文件
# NHANES官网链接：https://wwwn.cdc.gov/Nchs/Nhanes/2013-2014/BMX_H.htm

bmx.g <- read_xpt("2011-2012/Examination/bmx_g.xpt") #注意是 Examination 的类别
bmx.h <- read_xpt("2013-2014/Examination/bmx_h.xpt") #注意是 Examination 的类别

### 2.提取研究所需要的变量 BMI-BMXBMI; 腰围-BMXWAIST
bmx.data.file <- dplyr::bind_rows(list(bmx.g, bmx.h))
bmx.data <- bmx.data.file[,c('SEQN', 'BMXBMI', 'BMXWAIST')]

##### 1.5 饮食-L&Z 以及总卡路里数据提取(第1天 & 第2天) #####
# https://wwwn.cdc.gov/nchs/nhanes/2013-2014/DR1TOT_H.htm#WTDRD1
# 第1天
dr1tot.g <- read_xpt('2011-2012/Dietary/dr1tot_g.xpt')
dr1tot.h <- read_xpt('2013-2014/Dietary/dr1tot_h.xpt')
View(dr1tot.g)


dr1tot.data.file <- dplyr::bind_rows(list(dr1tot.g, dr1tot.h))
dr1tot.data <- dr1tot.data.file[,c('SEQN', 'DR1TLZ', 'DR1TKCAL')]


# 第2天
dr2tot.g <- read_xpt('2011-2012/Dietary/dr2tot_g.xpt')
dr2tot.h <- read_xpt('2013-2014/Dietary/dr2tot_h.xpt')
dr2tot.data.file <- dplyr::bind_rows(list(dr2tot.g, dr2tot.h))
dr2tot.data <- dr2tot.data.file[,c('SEQN', 'DR2TLZ', 'DR2TKCAL')]

# 合并第1天&第2天的变量
dr.data <- merge(dr2tot.data, dr1tot.data)
View(dr.data)
##### 1.6 膳食补充剂-L&Z数据提取(第1天 & 第2天) ######
# 第1天
ds1tot.g <- read_xpt('2011-2012/Dietary/ds1tot_g.xpt')
ds1tot.h <- read_xpt('2013-2014/Dietary/ds1tot_h.xpt')
# View(ds1tot.g)
ds1tot.data.file <- dplyr::bind_rows(list(ds1tot.g, ds1tot.h))
ds1tot.data <- ds1tot.data.file[,c('SEQN', 'DS1TLZ')]


# 第2天
ds2tot.g <- read_xpt('2011-2012/Dietary/ds2tot_g.xpt')
ds2tot.h <- read_xpt('2013-2014/Dietary/ds2tot_h.xpt')
ds2tot.data.file <- dplyr::bind_rows(list(ds2tot.g, ds2tot.h))
ds2tot.data <- ds2tot.data.file[,c('SEQN', 'DS2TLZ')]

# 合并第1天&第2天的变量
ds.data <- merge(ds1tot.data, ds2tot.data)
View(ds.data)
##### 1.7 CFQ-认知数据提取 #####
### 1.提取 Component文件
# NHANES官网链接：https://wwwn.cdc.gov/Nchs/Nhanes/2013-2014/CFQ_H.htm

cfq.g <- read_xpt('2011-2012/Questionnaire/cfq_g.xpt')
cfq.h <- read_xpt('2013-2014/Questionnaire/cfq_h.xpt')

### 2.提取研究所需要的变量
#CERAD 词汇第1次即刻记忆测验	CFDCST1
#CERAD 词汇第2次即刻记忆测验	CFDCST2
#CERAD 词汇第3次即刻记忆测验	CFDCST3
#CERAD 3次词汇即刻记忆测验总分	CERAD_1+CERAD_2+CERAD_3
#CERAD 延迟回忆评分	CFDCSR
#语言流畅性评分Animal Fluency	CFDAST
#数字符号替代测试(DSST)	CFDDS

cfq.data.file <- dplyr::bind_rows(list(cfq.g, cfq.h))
cfq.data <- cfq.data.file[,c('SEQN', 'CFDCST1', 'CFDCST2', 'CFDCST3',
                             'CFDCSR', 'CFDAST', 'CFDDS')]

#### 2 提取分析相关变量（权重等，暂时为复现paper结果而提取） ####
##### 2.1 权重变量 ##### 
# 找到权重变量
# DEMO->MEC-Dietary->MEC-CFQ
# 饮食中的权重：https://wwwn.cdc.gov/nchs/nhanes/2013-2014/DR1TOT_H.htm

# SMQ: https://wwwn.cdc.gov/Nchs/Nhanes/2013-2014/SMQ_H.htm
# SMQ 中没有权重信息：The NHANES full sample 2-Year MEC Exam Weights (WTMEC2YR) should be used to analyze the 2013-14 SMQ variables in conjunction with the laboratory measurements on tobacco exposure or other examination measurements.

# CFQ: https://wwwn.cdc.gov/Nchs/Nhanes/2011-2012/CFQ_G.htm
# CFQ 中没有权重信息：MEC exam sample weights should be used for the analyses 
# and are available on the NHANES website 

# dr1tot.data.add.weight <- dr1tot.data.file[,c('SEQN', 'WTDRD1', 'WTDR2D', 'DR1TLZ', 'DR1TKCAL')]
# dr2tot.data.add.weight <- dr2tot.data.file[,c('SEQN', 'WTDRD1', 'WTDR2D', 'DR2TLZ', 'DR2TKCAL')]
# dr.data.add.weight <- merge(dr1tot.data.add.weight, dr2tot.data.add.weight)

# View(dr.data.add.weight) #注意看数据中，weight = 0，没有对应的那一天的数据
# dim(dr.data.add.weight) #19151

weight.data <- dr1tot.data.file[,c('SEQN', 'WTDRD1')]
weight.data$WTDRD1 <- weight.data$WTDRD1/2 #权重的计算方式，需要按 Cycle 进行平均


##### 2.2 复杂抽样的其他变量-DEMO #####
survey.design.data <- demo.data.file[,c('SEQN', 'SDMVPSU', 'SDMVSTRA')]


#### 3.合并上述所有数据（把列拼接起来）-Output #### 
# 【4个# 用作第一层级】

##### 3.1 合并步骤1 & 2中提取的数据 #####
#demo.data, smq.data, alq.data, bmx.data, dr.data, ds.data, cfq.data, weight.data
output <- plyr::join_all(list(demo.data, smq.data, alq.data, bmx.data,
                        dr.data, ds.data, cfq.data, weight.data, survey.design.data),
                        by='SEQN', type='full')

# View(output) # 跑通到这一步，得到这篇文章所用到的全部原始数据, 19,931
dim(output) #[1] 19931    

##### 3.2 针对合并后的数据进行质控，确定人数能对得上 #####
# 1.根据 Table1 进行质控
data.age.60 <- subset.data.frame(output, RIDAGEYR >= 60 )
dim(data.age.60) #[1] 3632    

# 2. 确认下权重
which(!is.na(data.age.60$WTDRD1))

weight.non.na.index <- which(!is.na(data.age.60$WTDRD1))
sum(data.age.60$WTDRD1[weight.non.na.index]) # 计算权重-59659564 

#### 4.根据文章的筛选策略进行筛选 ####
# 原始筛选策略：The analyses presented here include a pooled sample of all individuals aged
# 60 years or older in each of the two survey cycles, with completed cognitive performance
# test results, dietary L and Z intake information, and information on important
# con-founders(age, sex, race/ethnicity, smoking, educational attainment).

# 注：这里不需要引号，e.g. 直接写 CFDCSR 就可以了（视频中存在错误）-12.06更新
paper.data <- subset.data.frame(output, RIDAGEYR >= 60 &
                                  (!is.na(CFDCSR)) & #CERAD 延迟回忆评分
                                  (!is.na(CFDAST)) & #语言流畅性评分Animal Fluency
                                  (!is.na(CFDDS)) & #数字符号替代测试(DSST)
                                  (!is.na(DR1TLZ))& #饮食-叶黄素 & 玉米黄质摄取
                                  (!is.na(DR1TKCAL))& #能量摄入总量-Energy (kcal)
                                  (!is.na(RIDAGEYR)) & #年龄
                                  (!is.na(RIAGENDR)) & #性别
                                  (!is.na(RIDRETH3)) & #种族
                                  (!is.na(DMDEDUC2)) & #教育程度
                                  (!is.na(SMQ020)) #是否吸烟至少100支
)
dim(paper.data)

# 注：这里不需要引号，e.g. 直接写 CFDCSR 就可以了（视频中存在错误）-12.06更新
paper.exclude.data <- subset.data.frame(output, RIDAGEYR >= 60 &(
                                          (is.na(CFDCSR)) | # CERAD 延迟回忆评分
                                          (is.na(CFDAST)) |# 语言流畅性评分Animal Fluency
                                          (is.na(CFDDS)) | # 数字符号替代测试(DSST)
                                          (is.na(DR1TLZ))| #饮食-叶黄素 & 玉米黄质摄取
                                          (is.na(DR1TKCAL))| #能力摄入总量-Energy (kcal)
                                          (is.na(RIDAGEYR)) | # 年龄
                                          (is.na(RIAGENDR)) | # 性别
                                          (is.na(RIDRETH3)) | # 种族
                                          (is.na(DMDEDUC2)) | # 教育程度
                                          (is.na(SMQ020))) # 是否吸烟至少100支
)
dim(paper.exclude.data)

index.DR.na <- (is.na(output$DR1TLZ) & is.na(output$DR2TLZ))
index.60.DR.na <- which(output$RIDAGEYR >=60 & index.DR.na)
length(index.60.DR.na) #[1] 564

index.CFDCSR.na <- is.na(output$CFDCSR)
index.60.DR.na <- which(output$RIDAGEYR >=60 & index.CFDCSR.na)
length(index.60.DR.na) #[1] 506


dim(paper.exclude.data) #[1] 916  
dim(paper.data) #[1] 2716   

# 将paper关注的2700+60岁以上有完整研究数据的人员，和排除的人员汇总在一起，并且打上标签
paper.exclude.data$Type <- 'exclude'
paper.data$Type <- 'include'
compare.data <- dplyr::bind_rows(paper.data, paper.exclude.data)



# tableby 函数使用的注意事项-1，分类变量可以在前面加上-factor；
# tableby 函数使用的注意事项-2，用于分队列的变量写在前面（Type变量）；
 
tab1 <- tableby(Type ~ factor(DMDEDUC2), data = compare.data) # 复现论文中的排除人群和纳入人群在人口学变量上的差异

summary(tab1)
View(paper.data)



